Source code for hysop.numerics.stencil.stencil_generator

# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.


"""Finite differences stencil inteface used in various numerical methods.

* :class:`~hysop.numerics.stencil.stencil.Stencil`
* :class:`~hysop.numerics.stencil.stencil_generator.StencilGenerator`

"""
import copy
import fractions
import gzip
import itertools as it
import math
import os

import numpy as np
import scipy as sp
import sympy as sm

try:
    import cPickle as pickle
except:
    import pickle

from hysop.numerics.stencil.stencil import CenteredStencil, Stencil
from hysop.tools.cache import load_data_from_cache, update_cache
from hysop.tools.io_utils import IO
from hysop.tools.misc import prod
from hysop.tools.numerics import F2Q, MPFR, MPQ, MPZ, mpq, mpqize, mpz
from hysop.tools.sympy_utils import (
    build_eqs_from_dicts,
    factor_split,
    tensor_symbol,
    tensor_xreplace,
)
from hysop.tools.htypes import extend_array

try:
    import flint

    has_flint = True
except ImportError:
    flint = None
    has_flint = False


[docs] class StencilGeneratorConfiguration: _dx = ("dx", "dy", "dz") def __init__(self): self.dim = 1 self.dtype = MPQ self.dx = sm.Symbol("dx") self.user_eqs = {} self.derivative = 1 self.order = 1 self.mask_type = StencilGenerator.CROSS self._mask = None self._format_inputs()
[docs] def copy(self): obj = StencilGeneratorConfiguration() for var in [ "dim", "dtype", "dx", "user_eqs", "derivative", "order", "mask_type", "_mask", ]: setattr(obj, var, copy.deepcopy(getattr(self, var))) return obj
[docs] def configure( self, dim=None, dtype=None, dx=None, user_eqs=None, derivative=None, order=None, mask=None, mask_type=None, ): """ Configure the stencil generator. """ self.dim = dim if (dim is not None) else self.dim self.dtype = dtype if (dtype is not None) else self.dtype self.dx = ( dx if (dx is not None) else (sm.Symbol(_) for _ in self._dx[: self.dim]) ) self.user_eqs = user_eqs if (user_eqs is not None) else self.user_eqs self.derivative = derivative if (derivative is not None) else self.derivative self.order = order if (order is not None) else self.order self.mask_type = mask_type if (mask_type is not None) else self.mask_type self._mask = mask if (mask is not None) else self._mask self._format_inputs() return self
[docs] def shape(self): if (self.derivative is None) or (self.order is None): raise RuntimeError("First set derivative and order with configure!") return self.derivative + self.order
[docs] def L(self, origin): StencilGeneratorConfiguration._check_origin(origin, self.dim) shape = self.shape() return origin
[docs] def R(self, origin): StencilGeneratorConfiguration._check_origin(origin, self.dim) shape = self.shape() return shape - origin - 1
[docs] def mask(self, origin, custom_mask=None): StencilGeneratorConfiguration._check_origin(origin, self.dim) self._check_and_raise() if self.mask == StencilGenerator.CUSTOM: if (custom_mask is None) and (self._mask is None): raise ValueError("Custom mask should not be None!") elif custom_mask is not None: mask = custom_mask.copy() self._mask = mask else: mask = self._mask else: mask = StencilGeneratorConfiguration._genmask( origin, self.shape(), self.mask_type ) return mask
[docs] def symbolic_derivatives(self, extra_size=0): df, df_vars = tensor_symbol(prefix="f", shape=self.shape() + extra_size) return df, df_vars
[docs] def symbolic_stencil(self, origin): StencilGeneratorConfiguration._check_origin(origin, self.dim) S, svars = tensor_symbol( prefix="S", mask=self.mask(origin), shape=self.shape(), origin=origin ) return S, svars
def _format_inputs(self): if ( isinstance(self.dx, sm.Symbol) or isinstance(self.dx, sm.Expr) or np.isscalar(self.dx) ): self.dx = np.asarray([self.dx] * self.dim) dim = self.dim self.derivative = extend_array(self.derivative, dim, dtype=int) self.order = extend_array(self.order, dim, dtype=int) self.dx = extend_array(self.dx, dim, dtype=object) def _check_and_raise(self): self._format_inputs() if self.dim < 1: raise ValueError("Dim < 1!") if self.dtype not in [MPQ, float, np.float32, np.float64]: raise TypeError(f"Invalid dtype {self.dtype}!") if not isinstance(self.dx, np.ndarray): raise TypeError("Invalid type for dx!") if self.dx.size != self.dim: raise ValueError("Invalid size for dx!") if not isinstance(self.user_eqs, dict): raise TypeError("Invalid type for user_eqs!") for k in self.user_eqs: if not isinstance(k, sm.Symbol): raise TypeError(f"Invalid type for key {k} in user_eqs!") if (self.derivative < 0).any(): raise ValueError("derivative < 0!") if (self.order < 1).any(): raise ValueError("order < 1!") @staticmethod def _check_origin(origin, dim): if not isinstance(origin, np.ndarray): raise TypeError("Origin is not a np.ndarray!") if origin.size != dim: raise ValueError("Bad dimensions for origin!") if (origin < 0).any(): raise ValueError("Origin component < 0!") @staticmethod def _genmask(origin, shape, mask): dim = origin.size if mask == StencilGenerator.CROSS: mask = np.zeros(shape, dtype=bool) for d in range(dim): access = [slice(origin[dd], origin[dd] + 1) for dd in range(dim)] access[d] = slice(None) access = tuple(access) mask[access] = True mask[tuple(access)] = True elif mask == StencilGenerator.DIAG: mask = np.ones(shape, dtype=bool) for d in range(dim): access = [slice(origin[dd], origin[dd] + 1) for dd in range(dim)] access[d] = slice(None) access = tuple(access) mask[access] = False mask[tuple(access)] = False mask[origin] = True elif mask == StencilGenerator.DENSE: mask = np.ones(shape, dtype=bool) else: raise NotImplementedError("Mask not implemented yet!") return mask def __str__(self): ss = """ StencilGeneratorConfiguration dim: {} dtype: {} dx: {} user_eqs: {} derivative: {} order: {} mask_type: {} mask: {} shape: {} """.format( self.dim, self.dtype, self.dx, self.user_eqs, self.derivative, self.order, self.mask_type, self._mask, self.shape(), ) return ss
[docs] class StencilGenerator: """ Generate various stencils. Results are cached and compressed locally. """ # Stencil masks DENSE = 0 CROSS = 1 DIAG = 2 CUSTOM = 99
[docs] @classmethod def cache_file(cls): _cache_dir = IO.cache_path() + "/numerics" _cache_file = _cache_dir + "/stencil.pklz" return _cache_file
def __init__(self, config=None, cache_override=False): """ StencilGenerator to generate stencils of dimension :attr:`dim` and type :attr:`dtype`. Results are cached and compressed locally. Attributes ---------- dim : int Dimension of the generated stencils. dtype : {np.float32, np.float64, :class:`MPQ`} Datatype of the generated stencils. Parameters ---------- dim : int Dimension of the generated stencils. dtype : {np.float32, np.float64, :class:`MPQ`} Datatype of the generated stencils. Other Parameters ---------------- cache_override : bool If set, overwrite already cached stencils. Raises ------ ValueError Raised if dim<1 or if dtype in [np.float32, np.float64] and dim!=1. TypeError Raised if dtype is not a valid type. Notes ----- Depending on data type :attr:`dtype` a different solver can be used. For fast generation of 1D stencils, use `np.float32` or `np.float64`. For n-dimentional exact fractional coefficients stencils use the default `mpq` type. Exact fractional results (mpq) are cached locally. Examples -------- Generate 1D approximative forward first derivative stencil of order 1 and 2. >>> import numpy as np >>> from hysop.numerics.stencil.stencil_generator import StencilGenerator >>> generator = StencilGenerator().configure(dim=1,dtype=np.float64) >>> s0 = generator.generate_approximative_stencil(origin=0, derivative=1, order=1) >>> s1 = generator.generate_approximative_stencil(origin=0, derivative=1, order=2) >>> print('{}\\n{}'.format(s0.coeffs,s1.coeffs)) [-1. 1.] [-1.5 2. -0.5] Generate 1D exact backward first derivative stencil of order 1 and 2. >>> from hysop.numerics.stencil.stencil_generator import StencilGenerator >>> generator = StencilGenerator().configure(dim=1) >>> s0 = generator.generate_exact_stencil(origin=-1, derivative=1, order=1) >>> s1 = generator.generate_exact_stencil(origin=-1, derivative=1, order=2) >>> print('{}\\n{}'.format(s0.coeffs,s1.coeffs)) [-1 1] [1/2 -2 3/2] Generate centered 2D Laplacian of order 2: >>> import sympy as sm >>> from hysop.numerics.stencil.stencil_generator import StencilGenerator >>> generator = StencilGenerator().configure(dim=2, derivative=2) >>> laplacian = generator.generate_exact_stencil(origin=1, order=2, dx=(sm.Symbol('dx'),)*2) >>> print(laplacian.coeffs) [[0 1 0] [1 -4 1] [0 1 0]] Similarly, with different custom dx: >>> laplacian = generator.generate_exact_stencil(origin=1, order=2) >>> print(laplacian.coeffs) [[0 dx**(-2) 0] [dy**(-2) (-2*dx**2 - 2*dy**2)/(dx**2*dy**2) dy**(-2)] [0 dx**(-2) 0]] See Also -------- :class:`Stencil` """ config = StencilGeneratorConfiguration() if (config is None) else config if not isinstance(config, StencilGeneratorConfiguration): raise TypeError( "config is not an instance of StencilGeneratorConfiguration!" ) self._config = config self._cache_override = cache_override @staticmethod def _hash_algorithm(): import hashlib return hashlib.new("sha512") @classmethod def _hash_equations(cls, eqs): hasher = cls._hash_algorithm() for e in eqs: hasher.update(str(e).encode("utf-8")) return hasher.hexdigest()
[docs] def get_config(self): """ Get current stencil generator configuration state. """ return self._config
[docs] def cache_override(self): """ Get caching override flag. """ return self._cache_override
[docs] def configure(self, config=None, **kargs): """ Push a new stencil generator configuration by keyword arguments (see StencilGeneratorConfiguration.configure() or by StencilGeneratorConfiguration instance. """ if config is not None: self._config = config self._config.configure(**kargs) return self
[docs] def generate_approximative_stencil(self, origin, **kargs): """ Generate a stencil using hardware floating point arithmetic (fast). dtype can be np.float16, np.float32, np.float64, MPQ """ config = self._config.copy() config.configure(**kargs) dim = config.dim dtype = config.dtype if dim != 1: raise ValueError("Bad dimension for approximation stencil generation!") if has_flint: solve_dtype = flint.fmpq elif dtype not in [np.float16, np.float32, np.float64]: solve_dtype = np.float64 else: solve_dtype = dtype dx = config.dx[0] k = config.derivative[0] order = config.order[0] N = config.shape()[0] origin = StencilGenerator._format_origin(origin, N) L = config.L(origin) R = config.R(origin) if k == 0: return Stencil([1], [0], 0, dx=dx, error=None) A = np.empty((N, N), dtype=solve_dtype) b = np.empty(N, dtype=solve_dtype) for i in range(N): b[i] = solve_dtype(int(i == k)) for j in range(N): A[i, j] = solve_dtype(int((j - origin) ** i)) try: if has_flint: coeffs = A.ravel() Afmpq = flint.fmpq_mat(*(A.shape + (coeffs,))) Afmpq_inv = Afmpq.inv() Ainv = np.asarray(Afmpq_inv.entries()).reshape(A.shape) S = Ainv.dot(b) else: S = sp.linalg.solve(A, b, overwrite_a=True, overwrite_b=True) except: print("\nError: Cannot generate stencil (singular system).\n") raise S *= math.factorial(k) actual_dtype = type(S.ravel()[0]) target_dtype = dtype if actual_dtype != target_dtype: if target_dtype in [np.float16, np.float32, np.float64]: if has_flint and (actual_dtype is flint.fmpq): def convert(x): return target_dtype(float(int(x.p)) / float(int(x.q))) S = np.vectorize(convert)(S) else: S = S.astype(target_dtype) elif target_dtype == MPQ: if has_flint and (actual_dtype is flint.fmpq): def convert(x): return mpq(mpz(x.p.str()), mpz(x.q.str())) else: def convert(x): frac = fractions.Fraction(x).limit_denominator((1 << 32) - 1) return mpq(frac.numerator, frac.denominator) S = np.vectorize(convert)(S) else: RuntimeError("Type conversion not implemented yet.") return Stencil(S, origin, order, factor=1 / (dx**k), dx=dx)
[docs] def generate_exact_stencil(self, origin, **kargs): """ Generate a stencil using quotients (slower). Yield exact solutions by using symbolic calculation and quotient arithmetic. Do not use on high order stencils. """ config = self._config.copy() config.configure(**kargs) config._check_and_raise() dx = config.dx dim = config.dim dtype = config.dtype order = config.order derivative = config.derivative N = config.shape() origin = StencilGenerator._format_origin(origin, N) L = config.L(origin) R = config.R(origin) df, df_vars = config.symbolic_derivatives(4) S, svars = config.symbolic_stencil(origin) user_eqs = config.user_eqs if all(derivative == 0): return Stencil([1], [0], 0, dx=dx, error=None) if len(user_eqs) == 0: for i, d in enumerate(derivative): if dim > 1: access = [slice(0, 1) for _ in range(dim)] access[i] = slice(d, d + 1) access = tuple(access) user_eqs[df[access].ravel()[0]] = 1 else: user_eqs[df[d]] = 1 def taylor(df, dx, N): expr = 0 for n in range(max(N)): expr += taylorn(df, dx, n, N) return expr def taylorn(df, dx, n, N): dim = dx.size def preficate(it): return (it <= N).all() and sum(it) == n nn = range(n + 1) itd = it.product(nn, repeat=dim) itd = filter(preficate, itd) expr = 0 for der in itd: expr += taylorn_term(df, dx, der) return expr def taylorn_term(df, dx, der): fact = np.asarray([math.factorial(d) for d in der]) return df[der] * prod((dx**der) / fact) expr = 0 for ids in np.ndindex(*N): offset = ids - origin expr += S[ids] * taylor(df, offset * dx, N) fact = factor_split(expr, df_vars) eqs = build_eqs_from_dicts(fact, user_eqs) svars = list([_ for _ in svars if not (_ == 0)]) stencil = Stencil(S, origin, order) i = 0 sol = {} cache_file = self.cache_file() while stencil.variables().intersection(svars): eqs_key = StencilGenerator._hash_equations(eqs) nsol = load_data_from_cache(cache_file, eqs_key) if (self._cache_override) or (nsol is None): nsol = sm.solve(eqs, svars) update_cache(cache_file, eqs_key, nsol) if not nsol: break sol.update(nsol) S = tensor_xreplace(S, sol) err = 0 while err == 0: for ids in np.ndindex(*N): offset = ids - origin err += S[ids] * taylorn(df, offset * dx, max(N) + i, N + i + 1) if err != 0: order += 1 i += 1 eqs = [] for k, eq in factor_split(err, df_vars).items(): if isinstance(eq, int): continue eq = sm.simplify(eq.xreplace(sol)) if eq.atoms(sm.Symbol).intersection(svars): eqs.append(eq) stencil = Stencil(S, origin, order, dx=dx, error=err) if (dx == dx[0]).all() and dx[0] != 1: stencil.refactor(dx[0] ** (-derivative[0])) return stencil
@staticmethod def _format_origin(origin, shape): dim = shape.size origin = ( np.asarray([origin] * dim, dtype=int) if np.isscalar(origin) else np.asarray(origin) ) if origin.size != dim: raise ValueError("Bad input dimensions!") return (origin + shape) % shape
[docs] class CenteredStencilGenerator(StencilGenerator): """ Generate various centered stencils. Results are cached and compressed locally. """ def __init__(self, config=None, cache_override=False): super().__init__(config, cache_override)
[docs] def generate_exact_stencil(self, **kargs): config = self._config.copy() config.configure(**kargs) shape = config.shape() origin = (shape - 1) // 2 stencil = super().generate_exact_stencil(origin=origin, **kargs) if stencil.is_centered(): return CenteredStencil.from_stencil(stencil) else: raise RuntimeError("Generated stencil is not centered!")
[docs] def generate_approximative_stencil(self, **kargs): config = self._config.copy() config.configure(**kargs) shape = config.shape() origin = (shape - 1) // 2 stencil = super().generate_approximative_stencil(origin, **kargs) if stencil.is_centered(): return CenteredStencil.from_stencil(stencil) else: raise RuntimeError(f"Generated stencil is not centered: {stencil.coeffs}")
if __name__ == "__main__": from hysop.tools.contexts import printoptions sg = StencilGenerator() with printoptions(precision=4): sg.configure(dim=1, derivative=0, order=2, dtype=np.float64) print("\ndim=1, 0th derivative, np.float64, approximative, shape=(5,):") for i in range(sg.get_config().shape()[0]): stencil = sg.generate_approximative_stencil(origin=i).reshape((5,)) print(f" origin: {i} => {stencil.factor} . {stencil.coeffs}") sg.configure(dim=1, derivative=1, order=2, dtype=np.float64) print('\ndim=1, 1st order first derivative, np.float64, approximative, shape=(5,):') for i in range(sg.get_config().shape()[0]): stencil = sg.generate_approximative_stencil(origin=i).reshape((5,)) print(f" origin: {i} => {stencil.factor} . {stencil.coeffs}") sg.configure(dim=1, derivative=2, order=2, dtype=np.float64) print('\ndim=1, 2nd order first derivative, np.float64, approximative:') for i in range(sg.get_config().shape()[0]): stencil = sg.generate_approximative_stencil(origin=i) print(f" origin: {i} => {stencil.factor} . {stencil.coeffs}") print("\ndim=1, 2nd order first derivative, exact:") sg.configure(dtype=MPQ, derivative=1) for i in range(sg.get_config().shape()[0]): stencil = sg.generate_exact_stencil(origin=i) print(f" origin: {i} => {stencil.factor} . {stencil.coeffs}") print("\ndim=1, 2nd order second derivative, exact:") sg.configure(dtype=MPQ, derivative=2) for i in range(sg.get_config().shape()[0]): stencil = sg.generate_exact_stencil(origin=i) print(f" origin: {i} => {stencil.factor} . {stencil.coeffs}") print("\n 2D Laplacian, 2nd order") sg.configure(dim=2, order=2) laplacian = sg.generate_exact_stencil(origin=1) print(laplacian.coeffs) laplacian = sg.generate_exact_stencil( origin=1, dx=[sm.Symbol("dy"), sm.Symbol("dx")] ) print("\n", laplacian.coeffs) df, _ = sg.get_config().symbolic_derivatives() user_eqs = {df[0][2]: sm.Symbol("Cx"), df[2][0]: sm.Symbol("Cy")} stencil = sg.generate_exact_stencil(origin=1, user_eqs=user_eqs) print("\n", stencil.coeffs) sg = CenteredStencilGenerator() sg.configure(derivative=2) print("\nCentered second order derivative stencils:") for i in range(1, 4): stencil = sg.generate_approximative_stencil(order=2*i, dtype=np.float16) print(f' {stencil.coeffs}') print() for i in range(1, 4): stencil = sg.generate_approximative_stencil(order=2 * i, dtype=MPQ) print(f" {stencil.coeffs}") print() for i in range(1, 4): stencil = sg.generate_exact_stencil(order=2*i) print(f' {stencil.coeffs}')